rm(list=ls())
library(tidyverse)
library(scales)
library(broom)
#df<-import data24.csv as df



# Figure 2

# Panel A ---- 
fit_c0 <- lm( 
  ftjudcourt~cneg+pneg+female+educ+pidc+pid_dk+age+black+hispanic,data=df)

panel_a <- tidy(fit_c0) |>
  filter(term %in% c("cneg","pneg")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "A. Court Evaluations",
    terms = ifelse(term == "pneg", "Police", "Courts")
  )

# Panel B ----
panel_b<-lm(ftjudcourt~cneg+pneg+pneg*pnc1+female+educ+pidc+pid_dk+age+black+hispanic,data=df)
panel_b1 <- tidy(panel_b) |>
  filter(term %in% c("pneg","cneg")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "B. Court Evaluations (conditional police effect)",
    terms = ifelse(term == "pneg", "Police", "Courts")
  )

b  <- coef(panel_b)
V  <- vcov(panel_b)

# names: "p_exp", "pnc", "p_exp:pnc"
Delta_hat <- b["pneg"] + b["pnc1"] + b["pneg:pnc1"]

Var_Delta <- V["pneg","pneg"] +
  V["pnc1","pnc1"] +
  V["pneg:pnc1","pneg:pnc1"] +
  2*( V["pneg","pnc1"] + V["pneg","pneg:pnc1"] + V["pnc1","pneg:pnc1"] )

SE_Delta  <- sqrt(Var_Delta)
CI        <- c(Delta_hat - 1.96*SE_Delta, Delta_hat + 1.96*SE_Delta)

#c(estimate = Delta_hat, se = SE_Delta, ci_low = CI[1], ci_high = CI[2])

md<-data.frame(term="pnc",estimate = Delta_hat, se = SE_Delta, ci_low = CI[1], ci_high = CI[2])
colnames(md)<-c("term","estimate","std.error","ci_low","ci_up")
md$Evaluations<-c("B. Court Evaluations (conditional police effect)");md$terms<-c("Police")

panel_b<-rbind(panel_b1,md) 



# Panel C ----

fit_p0 <- lm( 
  ftpolice~cneg+pneg+female+educ+pidc+pid_dk+age+black+hispanic,data=df)

panel_c <- tidy(fit_p0) |>
  filter(term %in% c("cneg","pneg")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "C. Police Evaluations",
    terms = ifelse(term == "pneg", "Police", "Courts")
  )






# Panel D ----
panel_d<-lm(ftpolice~cneg+pneg+cneg*cnp1+female+educ+pidc+pid_dk+age+black+hispanic,data=df)
panel_d1 <- tidy(panel_d) |>
  filter(term %in% c("pneg","cneg")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "D. Police Evaluations (conditional court effect)",
    terms = ifelse(term == "pneg", "Police", "Courts")
  )

b  <- coef(panel_d)
V  <- vcov(panel_d)

# names: "p_exp", "pnc", "p_exp:pnc"
Delta_hat <- b["cneg"] + b["cnp1"] + b["cneg:cnp1"]

Var_Delta <- V["cneg","cneg"] +
  V["cnp1","cnp1"] +
  V["cneg:cnp1","cneg:cnp1"] +
  2*( V["cneg","cnp1"] + V["cneg","cneg:cnp1"] + V["cnp1","cneg:cnp1"] )

SE_Delta  <- sqrt(Var_Delta)
CI        <- c(Delta_hat - 1.96*SE_Delta, Delta_hat + 1.96*SE_Delta)

#c(estimate = Delta_hat, se = SE_Delta, ci_low = CI[1], ci_high = CI[2])

md<-data.frame(term="cnp",estimate = Delta_hat, se = SE_Delta, ci_low = CI[1], ci_high = CI[2])
colnames(md)<-c("term","estimate","std.error","ci_low","ci_up")
md$Evaluations<-c("D. Police Evaluations (conditional court effect)");md$terms<-c("Courts")

panel_d<-rbind(panel_d1,md) 

# Plot figure 2 ----

plotdf <- dplyr::bind_rows(panel_a,panel_b,panel_c,panel_d) %>%
  dplyr::mutate(
    Evaluations = factor(
      Evaluations,
      levels = c(
        "A. Court Evaluations",
        "B. Court Evaluations (conditional police effect)",
        "D. Police Evaluations (conditional court effect)",
        "C. Police Evaluations"
      )
    ),
    terms = factor(terms, levels = c("Courts","Police"))
  )

mid_lab <- "B. Court Evaluations (conditional police effect)"
mid_lab1 <- "D. Police Evaluations (conditional court effect)"
plotdf <- plotdf %>%
  group_by(Evaluations, terms) %>%
  mutate(
    terms_new = case_when(
      Evaluations == mid_lab & terms == "Police" & row_number() == 1 ~ "Police",
      Evaluations == mid_lab & terms == "Police" & row_number() == 2 ~ "Police only",
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 1 ~ "Courts",
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 2 ~ "Courts only",
      TRUE ~ terms
    ),
    y_off = case_when(
      Evaluations == mid_lab & terms == "Police" & row_number() == 1 ~ -0.10,
      Evaluations == mid_lab & terms == "Police" & row_number() == 2 ~  0.10,
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 1 ~ -0.10,
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 2 ~  0.10,
      TRUE ~ 0
    )
  ) %>% ungroup()

plotdf$Evaluations<-factor(plotdf$Evaluations,levels=c("A. Court Evaluations","B. Court Evaluations (conditional police effect)","C. Police Evaluations","D. Police Evaluations (conditional court effect)"))

ggplot(
  plotdf,
  aes(y = terms, x = estimate, xmin = ci_low, xmax = ci_up,
      shape = terms_new, linetype = terms_new, color = terms_new)
) +
  geom_point(size = 3, position = position_nudge(y = plotdf$y_off)) +
  geom_errorbarh(height = .14, position = position_nudge(y = plotdf$y_off)) +
  geom_vline(xintercept = 0, linetype = 5, col = "grey20") +
  xlab("") + ylab("") +
  scale_x_continuous(limits = c(-16,3), n.breaks = 20, minor_breaks = 0) +
  scale_color_manual(
    values = c("Courts"="black","Police"="black","Police only"="black",
               "Police"="black","Police only"="black","Courts only"="black"),
    breaks = c("Courts","Police","Police only","Courts only"),   
    labels = c("Courts","Police","Police only","Courts only")
  ) +
  scale_shape_manual(
    values = c("Courts"=0,"Police"=1,"Police only"=17,"Courts only"=18),
    breaks = c("Courts","Police","Police only","Courts only"),
    labels = c("Courts","Police","Police only","Courts only")
  ) +
  scale_linetype_manual(
    values = c("Courts"="longdash","Police"="solid"
               ,"Police only"="dotted","Courts only"="dotted"),
    breaks = c("Courts","Police","Police only","Courts only"),
    labels = c("Courts","Police","Police only","Courts only")
  ) +
  guides(color = guide_legend(title = "Experiences"),
         shape = guide_legend(title = "Experiences"),
         linetype = guide_legend(title = "Experiences")) +
  facet_wrap(vars(Evaluations), nrow = 2, drop = FALSE) +
  theme_bw() +
  theme(legend.position = "bottom")



# Figure 3




# Figure 3

# Figure 3
# Panel A ---- 
fit_c0 <- lm( 
  ftjudcourt~cneg+pneg+neg_aff+nat_ft_avg+female+educ+pidc+pid_dk+age+black+hispanic,data=df)

panel_a <- tidy(fit_c0) |>
  filter(term %in% c("cneg","pneg")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "A. Court Evaluations",
    terms = ifelse(term == "pneg", "Police", "Courts")
  )

# Panel B ----
panel_b<-lm(ftjudcourt~cneg+pneg+pneg*pnc1+neg_aff+nat_ft_avg+female+educ+pidc+pid_dk+age+black+hispanic,data=df)
panel_b1 <- tidy(panel_b) |>
  filter(term %in% c("pneg","cneg")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "B. Court Evaluations (conditional police effect)",
    terms = ifelse(term == "pneg", "Police", "Courts")
  )

b  <- coef(panel_b)
V  <- vcov(panel_b)

# names: "p_exp", "pnc", "p_exp:pnc"
Delta_hat <- b["pneg"] + b["pnc1"] + b["pneg:pnc1"]

Var_Delta <- V["pneg","pneg"] +
  V["pnc1","pnc1"] +
  V["pneg:pnc1","pneg:pnc1"] +
  2*( V["pneg","pnc1"] + V["pneg","pneg:pnc1"] + V["pnc1","pneg:pnc1"] )

SE_Delta  <- sqrt(Var_Delta)
CI        <- c(Delta_hat - 1.96*SE_Delta, Delta_hat + 1.96*SE_Delta)

#c(estimate = Delta_hat, se = SE_Delta, ci_low = CI[1], ci_high = CI[2])

md<-data.frame(term="pnc",estimate = Delta_hat, se = SE_Delta, ci_low = CI[1], ci_high = CI[2])
colnames(md)<-c("term","estimate","std.error","ci_low","ci_up")
md$Evaluations<-c("B. Court Evaluations (conditional police effect)");md$terms<-c("Police")

panel_b<-rbind(panel_b1,md) 



# Panel C ----

fit_p0 <- lm( 
  ftpolice~cneg+pneg+neg_aff+nat_ft_avg+neg_aff+nat_ft_avg+female+educ+pidc+pid_dk+age+black+hispanic,data=df)

panel_c <- tidy(fit_p0) |>
  filter(term %in% c("cneg","pneg")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "C. Police Evaluations",
    terms = ifelse(term == "pneg", "Police", "Courts")
  )






# Panel D ----
panel_d<-lm(ftpolice~cneg+pneg+cneg*cnp1+neg_aff+nat_ft_avg+female+educ+pidc+pid_dk+age+black+hispanic,data=df)
panel_d1 <- tidy(panel_d) |>
  filter(term %in% c("pneg","cneg")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "D. Police Evaluations (conditional court effect)",
    terms = ifelse(term == "pneg", "Police", "Courts")
  )

b  <- coef(panel_d)
V  <- vcov(panel_d)

# names: "p_exp", "pnc", "p_exp:pnc"
Delta_hat <- b["cneg"] + b["cnp1"] + b["cneg:cnp1"]

Var_Delta <- V["cneg","cneg"] +
  V["cnp1","cnp1"] +
  V["cneg:cnp1","cneg:cnp1"] +
  2*( V["cneg","cnp1"] + V["cneg","cneg:cnp1"] + V["cnp1","cneg:cnp1"] )

SE_Delta  <- sqrt(Var_Delta)
CI        <- c(Delta_hat - 1.96*SE_Delta, Delta_hat + 1.96*SE_Delta)

#c(estimate = Delta_hat, se = SE_Delta, ci_low = CI[1], ci_high = CI[2])

md<-data.frame(term="cnp",estimate = Delta_hat, se = SE_Delta, ci_low = CI[1], ci_high = CI[2])
colnames(md)<-c("term","estimate","std.error","ci_low","ci_up")
md$Evaluations<-c("D. Police Evaluations (conditional court effect)");md$terms<-c("Courts")

panel_d<-rbind(panel_d1,md) 




# Plot figure 3 ----

plotdf <- dplyr::bind_rows(panel_a,panel_b,panel_c,panel_d) %>%
  dplyr::mutate(
    Evaluations = factor(
      Evaluations,
      levels = c(
        "A. Court Evaluations",
        "B. Court Evaluations (conditional police effect)",
        "D. Police Evaluations (conditional court effect)",
        "C. Police Evaluations"
      )
    ),
    terms = factor(terms, levels = c("Courts","Police"))
  )


mid_lab <- "B. Court Evaluations (conditional police effect)"
mid_lab1 <- "D. Police Evaluations (conditional court effect)"
plotdf <- plotdf %>%
  group_by(Evaluations, terms) %>%
  mutate(
    terms_new = case_when(
      Evaluations == mid_lab & terms == "Police" & row_number() == 1 ~ "Police",
      Evaluations == mid_lab & terms == "Police" & row_number() == 2 ~ "Police only",
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 1 ~ "Courts",
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 2 ~ "Courts only",
      TRUE ~ terms
    ),
    y_off = case_when(
      Evaluations == mid_lab & terms == "Police" & row_number() == 1 ~ -0.10,
      Evaluations == mid_lab & terms == "Police" & row_number() == 2 ~  0.10,
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 1 ~ -0.10,
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 2 ~  0.10,
      TRUE ~ 0
    )
  ) %>% ungroup()

plotdf$Evaluations<-factor(plotdf$Evaluations,levels=c("A. Court Evaluations","B. Court Evaluations (conditional police effect)","C. Police Evaluations","D. Police Evaluations (conditional court effect)"))

ggplot(
  plotdf,
  aes(y = terms, x = estimate, xmin = ci_low, xmax = ci_up,
      shape = terms_new, linetype = terms_new, color = terms_new)
) +
  geom_point(size = 3, position = position_nudge(y = plotdf$y_off)) +
  geom_errorbarh(height = .14, position = position_nudge(y = plotdf$y_off)) +
  geom_vline(xintercept = 0, linetype = 5, col = "grey20") +
  xlab("") + ylab("") +
  scale_x_continuous(limits = c(-16,3), n.breaks = 20, minor_breaks = 0) +
  scale_color_manual(
    values = c("Courts"="black","Police"="black","Police only"="black",
               "Police"="black","Police only"="black","Courts only"="black"),
    breaks = c("Courts","Police","Police only","Courts only"),   
    labels = c("Courts","Police","Police only","Courts only")
  ) +
  scale_shape_manual(
    values = c("Courts"=0,"Police"=1,"Police only"=17,"Courts only"=18),
    breaks = c("Courts","Police","Police only","Courts only"),
    labels = c("Courts","Police","Police only","Courts only")
  ) +
  scale_linetype_manual(
    values = c("Courts"="longdash","Police"="solid"
               ,"Police only"="dotted","Courts only"="dotted"),
    breaks = c("Courts","Police","Police only","Courts only"),
    labels = c("Courts","Police","Police only","Courts only")
  ) +
  guides(color = guide_legend(title = "Experiences"),
         shape = guide_legend(title = "Experiences"),
         linetype = guide_legend(title = "Experiences")) +
  facet_wrap(vars(Evaluations), nrow = 2, drop = FALSE) +
  theme_bw() +
  theme(legend.position = "bottom")



















